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Abstract 

We investigate the multi-dimensional Super Resolution problem on closed semi- 
algebraic domains for various sampling schemes such as Fourier or moments. We 
present a new semidefinite programming (SDP) formulation of the ti-minimization in 
the space of Radon measures in the multi-dimensional frame on semi-algebraic sets. 

While standard approaches have focused on SDP relaxations of the dual program 
(a popular approach is based on Gram matrix representations), this paper introduces 
an exact formulation of the primal ti-minimization exact recovery problem of Super 
Resolution that unleashes standard techniques (such as moment-sum-of-squares hier¬ 
archies) to overcome intrinsic limitations of previous works in the literature. Notably, 
we show that one can exactly solve the Super Resolution problem in dimension greater 
than 2 and for a large family of domains described by semi-algebraic sets. 

Keywords: super resolution; signed measure; semidefinite programming; total variation; 

semialgebraic domain. 


1 Introduction 


1.1 Super Resolution 

The early formulation of the Super Resolution problem can be identihed as the ability of 
faithfully reconstruct a high-dimensional sparse vector from the observation of a low-pass 
hlter. This situation models important applications in imaging spectroscopy [HGTB94], 
image processing [PPK03], radar imaging [OBP94], or astronomy [MM05]. As a theoretical 
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baseline, suppose that one wants to reconstruct a vector x* by solving a system of linear 
equations: 

Ax = b where x G , b := Ax* , b G and m N . (1) 

If the number s of non-vanishing components of x* is small and the matrix A enjoys some 
geometric property, namely the Null Space Property [CDD09], that depends only on its kernel 
and the sparsity s, then one can exactly reconstruct x* by minimizing the f'l-norm within the 
affine subspace of all solutions of the linear system. Conditions on m, N, s and properties 
have been extensively studied, see for example [Don06b], [CDD09] and reference therein. 
Less than ten year ago. Super Resolution has seeded the ideas of compressed sensing theory 
[Don06a], [CT06], [CT07]. In this theory, the matrix A is randomized and one is interested 
both in the construction of probability distributions allowing to show relevant properties, 
such as the Restricted Isometry Property [CT06], and the stability of the reconstruction 
process. This area of research is very fruitful and leads to many practical applications in 
signal and image processing, see for example [HS09], [TG07], or [CTL08]. To the best of our 
knowledge, the hrst mathematical works on Super Resolution are due to Donoho et al in the 
early ninety, see [DS89] and [Don92]. In these papers, the term Super Resolution appeared 
because the matrix A is related to a discretization of some Fourier transform. As a matter 
of fact, when inverting a discrete Fourier transform, the separation of two close spikes in a 
sparse signal is made possible by minimizing the £i-norm while a linear inversion method is 
not able to do so. Beside, in these years, many applied researchers were performing Fourier 
inversion of non negative sparse signals using some entropy regularization, see for example 
[SG87] and [Bur67]. At this time, the respective roles of sparsity and non negativity in 
regard of spikes separation from non linear Fourier inversion methods was not completely 
clear. An important step for understanding these roles has been taken by lifting the linear 
equation (1) up to the more abstract measure set up, see [Gas90], [GG96] and [DG96]: 

(aj,/r)= / ai{x)pi{dx) = bi where/r G ^(A)+and i = 1,..., m . (2) 

Jx 

Here, (aj(x)) is a vector of continuous function dehned on X, a given compact subset of M"', 
b = [bi) G M”* and is the set of all nonnegative Radon measures on X. In this frame, 

there exists very special points b* such that the set of all members of R?{X)^ satisfying (2) 
for b = b* reduces to a singleton {/ib*}- Furthermore, pb* is a discrete measure concentrated 
on very few points. Hence, if one deals with a b close to a point b* the set of all solutions 
of (2) (or (1) with positivity constraint and a design matrix A discretizing the vectorial 
function (aj(x))™;^) is very small. So that, two different methods for selecting a member 
of this set will lead to similar solutions. We refer to [Ana92], [GG96] and [DG96] where 
quantitative evaluations on the size of the set of solutions is performed and to [GG94] and 
[Lew96] for related evaluations in another context. Notice that the structure of points b* can 
be completely described in the case where the family of functions (oj) is a Ghebyshev system, 
T-system for short, this includes the case of discrete Fourier transform and moments, see 
[BE95] or [KS66] for an exhaustive overview on these systems of functions. A more involved 
situation is when (2) does not enjoy the non negativity assumption on the measure. This 
means that one wishes to solve the linear equation: 

(oj, P')= ai{x)pi{dx) = bi where /i G R?{X) and z = 1,..., m . (3) 

Jx 
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Here, ^{X) is the set of signed Radon measnres on X. This is the frame of the present paper. 
Snrprisingly, as shown by anthors of this paper [dCG12] and nnder some assnmptions on the 
family of fnnctions (oj), there exists pairs (b*,/ib*) £ IR™' x ^{X) such that b* = (aj,/ib*), 
for z = 1, ■ • ■ , m, and pb* is the unique solution of (3) minimizing the total variation norm 
with b = b*. Such /Xb* are sparse in the sense that they are measures with finite support. 

The study of the solutions to (1) and (3) uncovers that f'l-minimization faithfully recon¬ 
struct objects concentrated on very few points. However, the analysis in Super Resolution 
differs dramatically from Compressed Sensing. For instance, it is well known that sparse 
f'l-minimization cannot be successful in the ultrahigh-dimensional setting [Verl2] where 
N 3> exp(m). Indeed it was shown in [CGLP12] that one needs at least m > (cst)s log(A^/s) 
measures to faithfully uncover s sparse vectors by £i-minimization. Hence, the analysis of 
Compressed Sensing in terms of high-dimensional random geometry [CGLP12] cannot be 
extended to the space of measure. Moreover, observe that Compressed Sensing aims at re¬ 
covering a sparse signal from random projection while, in Super Resolution, the sampling 
scheme is deterministic. 

Admittedly their analysis differ but we can bridge the gap between (1) and (3) by considering 
their dual formulations. From the point of view of convex analysis, we see that the dual form 
of these linear programs aims at reconstructing a dual certihcate [CT06, dCG12], i.e. an 
£oo-constraint linear combination of (a*), where (oj) are the lines of A in (1) and a family 
of continuous functions in (3). As pointed out by authors of this paper [dCG12], a parallel 
between Compressed Sensing and Super Resolution exists where the lines of A are the 
evaluation of a vector of continuous function at some prescribed points. In this frame. Super 
Resolution can be seen as a Compressed Sensing problem where the dimension N goes to 
inhnity. This analogy persists with the notion of dual certificate, i.e. a solution to the 
following dual program (4). Indeed, the dual programs given by the constraints (1) and (3) 
(while minimizing the f'l-norm and the total variation norm respectively) share the same 
expression: 

sup b^u s.t. ||a'''u||oo< 1, (4) 

where a = A in Compressed Sensing (1) and a(a;) = {ai{x))^i is a vector of continuous 
function in Super Resolution (3). Define a dual certificate P as: 

P = a^u* , (5) 

where u* is a solution to the dual program (4). From the duality properties, note that P is 
a sub-gradient of the £i-norm at a solution to the primal program (3). Hence, we can ensure 
that a target measure /i* is a solution to (3) if we are able to construct a dual polynomial 
(5) that interpolates the phases of the weights of /z* at its support points. 

From a theoretical point of view, one of the main issue in Super Resolution consists in 
exhibiting such a dual certificate P. In the Fourier frame, notice that an important con¬ 
struction, for target discrete measures whose support satisfy a separable condition, is given 
in the fundamental paper [CFG14] where a huge step has been taken. Indeed, the authors 
are the first to give a sharp condition on the support points of the target measure in order to 
warrant the existence of a dual certificate. Moreover, their proof is based on interpolating, 
by a Jackson kernel, the phases of the weights of the target measure at its support points 
and, henceforth, explicitly construct a dual certihcate (5). In the Compressed Sensing frame. 
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observe that the same method has been investigated in [Kahll] nsing a Dirichlet kernel. In 
the present paper, we will not deal with this issne bnt rather with the practical resolntion 
of the convex program (4). 

In the Snper Resolution frame, remark that the program (4) has hnitely many variables but 
inhnitely many constraints. This last point can be a severe limitation in practice. As a 
matter of fact, a difficult task is to construct a tractable program that deals with the £oo- 
constraint of (4). Standard formulations [CFG14] are based on Gram matrix representations, 
see below. Incidentally, these procedures cannot be extended to dimension greater than 2 or 
to semi-algebraic domains X. To cope with this issue, we consider a new parametrization of 
the primal program based on works of the authors [LaslO]. Note that our method relies on 
inhnitely many parameters but relaxations involving only a hnite number of parameters are 
proved to lead to the exact solution of the primal program. 

1.2 Previous works 

During the last years, theoretical guarantees for exact recovery [BP13, dCG12, GFG14], 
bounds on the support recovery from inaccurate samplings [AdGG13, FG13], prediction 
of the Fourrier coefficient from noisy observations [TBR13], and noise robustness [DP13] 
have been showed. These works prove that discrete measures can be recovered, in a robust 
manner, from few samples using an £i-method. 

From a numerical point of view, a solution to £i-minimization is often computed using the 
dual program described by (4). Then, the ^oo-norm constraint of the dual program (4) is 
equivalently formulated as a nonnegative constraint on (trigonometric) polynomials. This 
point of view unleashes Gram matrix representations [GFG14] or Toeplitz matrix represen¬ 
tations [TBR13] to handle the constraint of non-negativity of (trigonometric) polynomials 
on domains. However, these formulations are limited to the frame of the real line and the 
torus in dimension one (see [Dum07] for instance) since they rely on the Fejer-Riesz theorem. 
As a matter of fact, the literature of Super Resolution has been focused on the fact that 
a semidefinite programming (SDP) formulation for the dual problem can be given as long 
as there exists a spectral factorization for a globally nonnegative trigonometric polynomial 
[GFG14] (Bounded Real Lemma using Fejer-Riesz theorem) or a spectral decomposition for 
semi-dehnite Toeplitz matrices [TBR13] (Garatheodory-Toeplitz theorem). Hence, except 
on the real line and the torus in dimension one, there is no exact SDP formulation of the 
Super Resolution problem for the dual form. However, relaxed SDP versions of the dual 
form in dimension greater than 2 are discussed in [XGV-l-13] and they have been used on 
the 2-sphere in [BDF15, BDF14]. 

1.3 Contribution 

To the best of our knowledge, the present paper is the first to overcome this limitation 
and expand the scope of Super Resolution implementation to the multi-dimensional frame 
in general basic semi-algebraic domains. Indeed, we operate a smart method to tackle 
numerically the solution of equation (3) with minimal £i-norm. This method uses both a 
re-parametrization in terms of moment sequences [LaslO] and the so-called sum-of-squares 
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(SOS) decompositions of nonnegative multidimensional polynomials, used widely in systems 
control theory during the last decade, see for example [HG05]. Notice that, in the scope of 
Super Resolution, this technique is new and, contrary to other approaches, focuses on the 
primal program through a truncation of the moment sequences. 

More precisely, given the real numbers 6*, i = 1,... ,m, consider the inhnite-dimensional 
optimization problem: 

inf II/xIItv 

s.t. {ai,fi) = bi, i = l,...,m (6) 

H G ^{X), 

where || . \\tv is the total variation norm of measures (to be dehned later). Notice that, under 
standard assumptions, problem (6) is feasible. That is 

3/r e ^{X) such that for z = 1,..., m , (oj, /i) = 6j. (7) 

Our main contribution concerns the numerical resolution of the total variation minimization 
problem (6). We extend the univariate (n = 1) trigonometric SDP formulation of [CFG14] 
to a much more general SDP formulation in dimension n > 2, for measures supported on 
basic semialgebraic sets. 

To this end, we use the Jordan decomposition of the signed measure /x = /i+ — as a 
difference of two nonnegative measures supported on X and we follow [LaslO] to dehne a 
hierarchy of hnite-dimensional primal-dual SDP problems: 

• the primal problems correspond to SDP relaxations of the conditions that must satisfy 
hnitely many moments of the two nonnegative measures on X; 

• the dual problems correspond to SDP strengthenings using SOS multipliers of the 
conditions that two distinguished polynomials are nonnegative on X. 

The moment-SOS hierarchy is indexed by an integer k, called relaxation order, which is the 
(half of the) number of moments used to represent the measures in the primal problem, or 
equivalently, the (half of the) degree of the SOS representations of the polynomials in the 
dual problem. The larger is the relaxation order k, the larger is the size of the SDP problems, 
the number of variables and constraints growing polynomially in 0{kP). 

The primal SDP problem features the matrices of moments of the two nonnegative measures. 
If the rank of each moment matrix, as a function of fc, stabilizes to a certain constant value, 
then the corresponding measure is atomic, with the number of atoms equal to the rank. 
Therefore, the total variation minimization problem (6) has been solved succesfully, and this 
is certihed by the polynomials solving the SDP problem. Numerical linear algebra can then 
be used to retrieve the support of the optimal measure. 

In the sequel, we present some examples for which our method is the hrst to give an SDP 
formulation of the Super Resolution phenomena. As a matter of fact, our procedure en¬ 
compasses a larger class of measurements than the class of standard moments discussed 
previously. Our numerical experiments are carried out with the Matlab interface GloptiPoly 
3 which is designed to generate semidehnite relaxations of measure LP problems with poly¬ 
nomial data. So we assume that the functions ai{x) in LP problem (8) are multivariate 
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1 Ti 

polynomials, and for notational simplicity, we let ai(x) := a;"* = ■ ■ ■ Xn where a, G N” 

are given. Note that the choice of monomials is only motivated for notational simplicity, 
and that other choices of polynomials (e.g. Chebyshev polynomials) are typically preferable 
nnmerically^. SDP relaxations are then solved with SeDnMi or MOSEK, implementations 
of a primal-dnal interior-point algorithm. For reprodncibility pnrposes, onr Matlab codes 
(nsing the pnblic-domain interface GloptiPoly and the SDP solver SeDnMi) of the nnmerical 
examples presented next are available for download at 

homepages.laas.fr/henrion/software/tvsdp.tar.gz 


1.3.1 Disconnected domain 



Fignre 1: Degree 9 polynomial certihcate for the nnivariate example, with 2 points (red) in 
the snpport of the positive part, and 1 point (bine) in the snpport of the negative part of 
the optimal measnre. 

We want to recover the measnre: 

h ■= <^-3/4 + <^1/2 ~ <^1/8 

on the disconnected set X := [—1,—1/2] U [0,1] which can be modeled as the polynomial 
snperlevel set X = {r G M : gi{x) > 0} for the choice: 

gi{x) := —{x -I- l)(r -I- l/2)x{x — 1). 

numerical analysis of the impact of the basis is however out of the scope of our work. 


6 















In LP (8) we let ai{x) := x* and bi := (—3/4)* + (1/2)* — (1/8)* for i = 0,1, 2 ..., 9. Solving 
the SDP relaxation of order /c = 5 on our standard PC takes 0.2 seconds, and optimality 
is certified from the solution of the primal moment problem with a rank 2 moment matrix 
for /r+ and a rank 1 moment matrix for /i_, from which the 3 points can be extracted 
using numerical linear algebra. On Figure 1 we represent the degree 9 polynomial 
certifying optimality, constructed from the solution of the dual SOS problem. Indeed we can 
check that the polynomial attains the value +1 at the points x = —3/4, and x = 1/2 (in 
red), it attains the value —1 at the point x = 1/4 (in blue), while taking values between — 1 
and +1 on X. Notice in particular that the polynomial is larger than +1 around x = —1/4, 
but this point is not in X. 

1.3.2 Low-pass filters in dimension greater than 3 

In the Fourier frame, the recent SDP formulations of £i-minimization in the space of complex 
valued measures are based on the Fejer-Riesz theorem. As a consequence, they cannot handle 
dimensions greater than 3. Observe that our procedure can bypass this limitation. For sake 
of readability, we present an example in dimension 2 although it can be extended to any 
dimension. We want to recover the measure 

h •= ^(-1/2,1/2) + <^(1/2,-1/2) + <^{1/2,1/2) + <5(0,0) ~ <5(0,-1/2) — <5(1/2,0) 

on the box X := [—1,1]^, from the knowledge of moments of degree up to 12, i.e. aj(x) = x* 
for z = 0,1,..., 12. Solving the SDP relaxation of order fc = 6 on our standard PC takes 
less than 3 seconds, and optimality is certihed from the solution of the primal moment 
problem with a rank 4 moment matrix for /r_|_ for and a rank 2 moment matrix for /r_ 
from which the 6 points of the support of the optimal measure /i can be extracted using 
numerical linear algebra with a relative accuracy around 10“®. On Figure 2 we represent 
the degree 12 polynomial certifying optimality, constructed from the solution of the dual 
SOS problem. Indeed we can check that the polynomial attains the value +1 at the 3 
points X G {(—1/2,1/2), (1/2,—1/2), (0, 0)}, it attains the value —1 at the 2 points x G 
{(0, —1/2), (1/2,0)} (in blue), while taking values between —1 and +1 on X. 

1.3.3 Localization of points on the sphere 

Recent extensions of Super Resolution to spike deconvolution on the 2-sphere from spherical 
harmonic measurements has been investigated in [BDF15, BDF14]. In these paper, the au¬ 
thors give a sufficient condition for exact recovery using f'l-minimisation and they investigate 
spikes localization when the measurements are perturbed by additive noise. 

From a numerical point of view, they used a relaxed version of the dual program (bounded 
real lemma in dimension d = 3 and a Gram representation of the £oo-constraint appearing 
in the dual). Our work naturally extends to this frame and provides an exact formulation 
of the primal form. 

For sake of numerical code simplicity, we have considered polynomials on restricted to 
the domain X given by the 2-sphere (note that one could have used homogenous spherical 
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Figure 2: Degree 12 polynomial certificate for the bivariate example, with 4 points (red) in 
the support of the positive part, and 2 points (blue) in the support of the negative part of 
the optimal measure. 


harmonics instead as in [BDF15]). We want to recover the measure 


/i : = 


<^(1,0,0) + <^(0,1,0) + <5(0,0,1) ~ <5^y2^y2 Qj — q yl) 



y? ys-, 
2 ’ 2 


that is supported on the positive orthant just for better visualization purposes. In LP (8), the 
Oj consist of 3-variate monomials of degree up to 5, i.e. m = 56. Solving the SDP relaxation 
of order k =6 on our standard PC takes less than 20 seconds, and optimality is certihed 
from the solution of the primal moment problem with a rank 3 moment matrix for /i_|_ and 
a rank 3 moment matrix for /r_, from which the 6 points can be extracted using numerical 
linear algebra. On Figure 3 we represent the degree 6 polynomial certifying optimality on 
the 2-sphere, constructed from the solution of the dual SOS problem. Indeed we can check 
that the polynomial attains the value -|-1 at the 3 prescribed points, it attains the value — 1 
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Figure 3: Degree 6 polynomial certificate for the sphere example, with 3 points (red) in the 
support of the positive part, and 3 points (blue) in the support of the negative part of the 
optimal measure. 



(in red), at the 3 others prescribed points (in blue), while taking values between —1 and +1 
on X. 

From a theoretical point of view, the minimal separation condition appearing in [BDF15] 
requires a degree 2 polynomial. So our example satisfy the sufficient aforementioned condi¬ 
tion. 


2 Primal and dual LP formulation 

2.1 General model and notation 

Let n be a positive integer. Denote by ]R[a:] the set of all polynomials on M”, and for d G N, 
the set of all polynomials on M"" with degree not greater than d. Further, we use the 
following notation: 
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• X C M”, is a given closed basic semi-algebraic set: 


X := {x : gj{x) > 0, j = 1,..., nx} 


where G M[a;], j = 1,... ,nx, are given polynomials whose degrees are denoted by 
dj, j = 1,... ,nx- It is assumed that X is compact with an algebraic certihcate of 
compactness. For example, one of the polynomial inequalities gj{x) > 0 should be of 
the form; 

n 

R'^ >0, 

i=l 

for R a sufficiently large constant. Let gx ■= {gj)j=i,...,nx- 

• Let a = be a linearly independent family of polynomials of degree at most d on 

X. Notice that m < {1 + d)'^. 

• For monomials we use the multi-index notation 


n 



i=i 


for every x = (xi, ..., Xn) G M” and a = (ai, ..., an) G N”. 


• ^{X), the space of continuous functions on X, a Banach space when equipped with 
the sup-norm: 


ll/ll = sup |/(a;)|. 

x&X 


• ^{X), the space of signed Radon measures on X, a Banach space isometrically iso¬ 
morphic to the topological dual '^(X)* when equipped with the total variation norm; 


llhllrv = lh(^)l , 

' Eev 


where the supremum is taken over all partitions V oi X into a hnite number of disjoint 
measurable subsets. 


• ^(X)+ C ^(X) and ^(X)+ C ^(X) the respective positive cones of nonnegative 
continuous functions on X and (nonnegative) Radon measures on X. We use the 
standard notation / > 0 and /i > 0 for membership in ^(X)_|_ and <^^(X)_|_, respectively. 


• To denote the integration of a function against a measure, we use the duality bracket: 


if, h) 



for all / G ^(X), and g, G ^(X). 
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With the usual Jordan decomposition: 


/i = /i+ - /r- 

into a sum of two nonnegative Borel measures /x+,/i_, the optimization problem (6) can be 
rewritten equivalently as a linear programming (LP) problem in the convex cone + , 

namely: 

p* = inf 

s.t. {tti, p+) - {tti, fi_) = bi, i = 

/i+ e ^(X)+ ^ ^ 

/i_ e ^(X)+. 

If a = and b = ^ problem (8) is the dual of the following 

LP problem in the convex cone ^(X)+: 

d* = sup b^u 

s.t. Zj^{x) ■.= 1 + sJ{x)\i (9) 

zJiyX) := 1 — a'''(a;)u G '^(X)+ 

where the maximization is w.r.t. u = G Remark that LP problem (9) can 

be also written as: 

d* = sup b^u 

s.t. ||a^(a;)u||oo < 1- 

Lemma 1 There is no duality gap between primal LP (8) and dual LP (9), i.e. p* = d*. 
Proof: Define the vector r{fi^,fi_) G by: 

r(/i_|_, p_ ) . ((1) p+) p T—) y (®1) T+) (®1) T—)y ■ ■ ■ {pmy T+) i^my L— )) 

and the set 

R := {r(/r+,/i_) : (/i+, /i_) G ^(X)+ X ^(X)+} C 

By [Bar02, Theorem 7.2], p* = d* provided that p* is finite and R is closed. Finiteness of p* 
follows from Assumption 7 and nonnegativity of the objective function (l,/i+) + (l,/i_). To 
prove that R is closed we have to show that for any sequence (/i”, ^ x.^(X)+ 

such that r(/i”,/i”) —)■ s G as n —)■ oo, one has s = r(/x, z/) for some finite measures 

li,iz E t^(X)+. Since the supports of all the measures are contained in a compact set, 
and since (l,/i+) + (1,/iL) —t so all measures p\y pLi are uniformly bounded. Therefore, 
from the weak-* compactness (and weak-* sequential compactness) of the unit ball (Banach- 
Alaoglu’s Theorem), there is a subsequence p^)km that converges weakly-* to an 

element (/r, z/) G ^(X)+ x ^(X)+. In particular, as all a* are continuous, 

lim r{p^, p^) = r(/r, z/), 

k—^oo 


which proves that R is closed. □ 

Lemma 2 For the dual LP problem (9) the supremum is attained. 
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Proof: The feasibility set 


f/ := {u e : ||a’^(x)u||oo < 1} 

of the LP problem (9) is a closed convex subset of a finite-dimensional Euclidean space, 
and it contains the origin. Since the objective function in LP (9) is continuous on U, the 
optimum is attained if U is bounded. Suppose that U is not bounded. Then there exists a 
sequence (u„)ngN C U such that ||u„|| —)■ oo as n — )■ oo. Write = A„v„, with ||v„|| = 1. 
Notice that 0 < A„ —)• oo and v„ G 17 because 0 E U and U is convex. Then 

||a ||a (3:)A,2^nIIoo Ajj||a ( 3 :)v^||qo ^ 1 ; 

so that ||a'''(a;)v„||oo < A^^ —)■ 0 as n —)■ cx). Since ||v„|| = 1, there exists a subsequence 
Uk and V with ||v|| = 1 such that —)■ v as /c —)■ cxd and ||a'''(a;)v||oo = 0. By linear 
independence of a, this implies that v = 0 , a contradiction. □ 

As a consequence of strong duality of Lemma 1, for any optimal primal-dual pair (/i, m) we 
have the complementarity conditions 

(^+,/i+) = 0, (^_,/i_) = 0 

implying jointly with Lemma 2 that 

spt C {x E X : a"'" (3:) u = 1 } 


and 


spt /i_ C {x E X : a'''(x) u = —1} 

for some continuous function x h->■ a^(3:)u, and where spt/x denotes the support of /x, that 
is, the smallest closed set S' C M" such that /i(]R" \ S') = 0. 


Lemma 3 Problem (6) has an optimal atomic measure supported on at most 2{m+l) points. 

Proof: Let /i+ be a nonnegative measure solving problem (8) and let 6 +j := (aj,/x+), with 
Oo := 1, X = 0,1,..., m. If 6 + = 0 then /i+ = 0 is trivially atomic (with no atoms), so 
assume 6+0 7 ^ 0 , and consider the probability measure /x+ := /i+/ 6 +o which satisfies the m 
equality constraints {ai,fi+) = 6+i := 6 +j/ 6 +o, i = I, ■ ■ ■ ,m. From [Bar02, Proposition 9.4] 
there exists a probability measure /x+ satisfying the same equality constraints {at, fi+) = 6 +j 
and which is supported on (at most) m -|- 1 points of X. The same reasoning can be applied 
to any nonnegative measure /i_ solving problem (8), which has a discrete counterpart /x_ 
supported on (at most) m -|- 1 points of X. The result follows by considering the union of 
these two discrete supports, which consists of (at most) 2(m -|- 1) points of X. □ 


3 Primal and dual SDP formulation 

Problem (8) is an instance of a generalized moment problem. As such it can be solved by 
a converging hierarchy of finite-dimensional primal-dual semidefinite programming (SDP) 
problems, as described comprehensively in [LaslOj. In the sequel, we extract the key instru¬ 
mental ingredients to the construction of the hierarchy. 
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3.1 Primal moment SDP 


Recall from paragraph 2.1 that 

X := {x e M'" : gj{x) > 0, j = 1 ,... ,nx} 

is a basic semi-algebraic set with a straighforward certihcate of compactness, and let gx '■= 
{gj)j=i,...,nx denote its dehning polynomials. Given a measure /i G ^+(X), the real number 

?/„:=(x",/i) (10) 

is called its moment of order a G N”. Conversely, given a real valued sequence y := (|/a)aeN’», 
if identity (10) holds for all a G N"’, we say that y has a representing measure y G 
Equivalently, sequence y belongs to the inhnite-dimensional moment cone 

^(X) := {(ya)aem ■ l/a = (a;“,h), /i G ^(X)+}. 

In the sequel we describe a procedure to approximate this convex cone. 

Given k G N, let M[x]fc denote the space of real polynomials of degree at most k. Let 
us identify a polynomial p{x) = ^ with its vector p of coefficients in the 

monomial basis. Dehne the Riesz functional iy as the linear functional acting on polynomials 
as follows: p G M[x]a: e-)■ iy{p) = J^aPaVa = P^P ^ 1^- Notc that if sequence y has a 
representing measure /i, then fy(p) = {p,p)- Dehne the moment matrix of order k as the 
Gram matrix of the quadratic form p G h-)■ iy{p‘^) G M, i.e. the matrix Mk{y) such 

that iy{p^) = p~^Mk{y)p. By construction this matrix is symmetric and linear in y. Given 
a polynomial g G M[x], dehne its localizing matrix of order k as the Gram matrix of the 
quadratic form p G M.[x]k e->■ iy{gp‘^) G M, i.e. the matrix Mk{gy) such that ly{gp^) = 
p~^Mk{gy)p. By construction this matrix is symmetric and linear in y. For j = 1,, nx, let 
kj denote the smallest integer not less than half the degree of polynomial gj, and let kx '■= 
ma.x{l, ki,... ,kn^}. With these notations, and for k > kx, dehne the hnite-dimensional 
moment cone 

.y^kigx) ■= {{ya)\a\< 2 k ■ Mk{y) h 0, Mk-kj{gjy) ^ 0, j = 1, . . .,nx} 

where ^ 0 means positive semidehnite. 

Let y+ resp. ?/_ denote the sequence of moments 

y+a ■= j x“/i+(dx), y-a ■= j x°‘iJ,-{dx) 

of /i+ (resp. /i_), indexed by a G N”. Primal measure LP (8) can be written as a primal 
moment LP; 

p* = min y+o + y-o 

s-t. A{y+,y_) = h 

y+ e .^(X) 

y_ G .^(X) 
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where the linear system of equations K{y^y-) = b models the linear moment constraints. 
The moment relaxation of order k > maxj/cx, d} of the primal moment LP then reads: 

pI = min y+o + y_o 

s-t. A{y+,y_)=h . . 

y+ e ^fc(gx) 
y- G ^fc(gx) 

where the minimization is w.r.t. a vector (|/+,|/_) of moments of degree at most 2k. For 
hxed k, problem (11) is a hnite-dimensional linear programming problem in the convex 
cone of positive semidehnite matrices, i.e. an SDP problem. When k varies, the number 
of moments, as well as the size of the moment and localizing matrices in problem (11) are 
binomial coefficients growing in 0{k^). 

It can be shown that {pD is a monotonically nondecreasing converging sequence of lower 
bounds on p*, i.e. pl_^_i > pi and lim^^ooPfe = P* ■ However, in the context of solving LP (8), 
a more relevant result is the following: 

Theorem 1 For a given relaxation order k > maxjd, fcx}? l^t {yXiU*-) denote the solution 
of the moment SDP (11). If 


rank Mk-kx ivl) = rank Mk{yX) and rank Mk-kx iV-) = rank Mk{y*_) (12) 

then pI = p* and LP (8) has an optimal solution with pf (resp. p*_) atomic 

supported at r+ := rankMfc(yi^), (resp. r_ := rsuak Mk{yl.)) points. 

Proof: By [LaslO, Theorem 3.11], yl_ (resp. yf) is the vector of moments up to order 2k, 
of a measure (resp. p*_) supported on rank Mk{y\) points (resp. rank Mk{y*_)) points 
of X. Therefore p*_) is a feasible solution of (8) with value pk < p*, which proves that 
p*-) is an optimal solution of (8) and pk = p*. □ 

Given a moment matrix Mk{yl.) satisfying the rank constraint of Theorem 1, there is a 
numerical linear algebra algorithm that extracts the r+ points of the support of the cor¬ 
responding atomic measure p\, and similarly for p*_. The algorithm is described e.g. in 
[LaslO, Section 4.3] and it is implemented in the Matlab toolbox GloptiPoly 3. 

A certificate of optimality can be obtained by solving the dual problem to primal SDP 
problem (11), and this is described next. 


3.2 Dual SOS SDP 

For a given integer k, let Ii[x]k C M[a;] 2 A: denote the space of SOS (sums of squares) polyno¬ 
mials of degree at most 2k. If p G E[a;]fc this means that there exists qj G M[x]fc, j = 1,... ,np, 
such that p = Qj- Let 

:= {p G M[x] : p{x) > 0, Vx G X} 
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denote the infinite-dimensional cone of nonnegative polynomials on X, and for k > kx, 
define the finite-dimensional SOS cone, also called quadratic module 

nx 

^k{gx) ■= {Po + Pj e S[a;]fc_fc., j = 0,1,.. .,nx} C R[x]2k- 

Under the above assumptions on the polynomial family gx defining X, from Putinar’s the¬ 
orem, see e.g. [LaslO, Theorems 2.14 and 3.8], it holds that ^{X) fl M[a;]K is the closure 
of {Uk>kx'^k{Sx)) n M[a;]K for all n > kx- Observe also that ^{X) (resp. ^^(gx)) is the 
dual cone to ^(X) (resp. ^^(gx))- Whereas testing whether a given polynomial belongs 
to ^(X) is a difficult task, testing whether a given polynomial belongs to for a 

fixed k, amounts to solving an SDP problem. 

Dual continuous function LP (9) can be written as a positive polynomial LP 

d* = max b^u 

s.t. 1-I-a'''(a;) u G ^(X) 

1 — a'''(a;) u e ^(X) 

and its SOS strengthening of order k > kx reads: 

dl = max b^u 

s.t. 1-ha'^(a;)u e ^fc(gx) (13) 

1 -a’^(a;)u e ^k{gx) 

where the maximization is w.r.t. a vector u G M™. It turns out that this is SOS problem 
(13) is an SDP problem dual to the moment problem (11): 

Lemma 4 There is no duality gap between SDP problems (11) and (13), i.e. = dl, and 
both (11) and (13) have an optimal solution. 

Proof: We first show that (11) has an optimal solution. Recall that one of constraints 
gj{x) > 0 that define X states that M — ||a;|p > 0 for some M > 1. From the constraint 
Mk-kj{gjy+) P 0 one deduces that iy+{M — xj) > 0, and iyj^{Mxl — x\^‘^) > 0 for every 
t = 1,... ,2k — 2. Hence £y+{xj'^) < for every i = 1,... ,n. With similar arguments, 

Ey-^xf^) < M^y^o for every i = 1,... ,n. By [LaslO, Proposition 3.6] \y+a\ < M^y+o and 
\y-a\ < for all a G Next, in a minimizing sequence (p+,pi), s G N, of (11) 

one has -|- pig < y^^ + pig =: p for all s, and so \y^a\ — M^P and < M^p for 

all a G N 2 fc, and all s = 1,.... From this we deduce that there is a subsequence (|/+,|/!.‘), 
t eN, that converges to some (p^, yL) as t —)■ oo, with value y^Q + y*_Q = pi- In addition by 
a simple continuity argument, Mk{y\) P 0 and Mk-kj{gj 2/+) ^ 0, j = 1,... ,nx. Similarly 
Mk{y*_) P 0 and Mk-kj{,gjy*-) ^ 0, j = 1,... ,nx, which proves that {y\,y*_) is an optimal 
solution of (11). 

Next, the set of optimal solutions y* := {(?/+,pi)} of (11) is compact. This follows from 
\y*+a\ ^ ^^y +0 ^ M’^Pk |pl„| < M^y*_Q < M^pI for all a G And so every sequence 
in y* has a converging subsequence. From [Bar02, Chapter IV. Theorem 7.2] one also deduces 
that there is no duality gap between (11) and (13). 
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It remains to prove that (13) has an optimal solution. Consider a maximizing sequence 
(ui)tgN, with b^Ui pI = p* as t ^ oo. By feasibility in (13), one has ||a(a;)^ut||oo < 1 
for all t and therefore (u^) C W := {u e M"' : ||a(x)^u||oo < 1} and U is compact (see 
the proof of Lemma 2). Therefore there exists u* E U and a subsequence such that 

— )■ u* G W as £ —)■ oo. In particular b^u* = p*. Moreover, since by Lemma 8 in the 
Appendix the convex cone ^k{Sx) is closed, 1 — a(x)^Ut^ —)■ 1 — a(x)^u* G ^k{Sx), which 
proves that u* is an optimal solution of (13). □ 

Assume that the rank conditions of Theorem 1 is satished at some relaxation order k, and let 
P-) denote the atomic measures optimal for problem (8), obtained from the solution of 
the primal SDP problem (11). Let u* denote an optimal solution of the dual SDP problem 
(13). The duality result of Lemma 4 implies that 

Supp /i^ C {x G X : a'''(x)u* = 1} 

and 

Supp p*_ C {x E X a''^(x)u* = —1} 

so that the polynomial a^(x)u* can be used as a certihcate of optimality. We formulate this 
in the following dual to Theorem 1. 

Lemma 5 Assume that the rank conditions (12) of Theorem 1 hold. Let us denote by u* 
the optimal solution of SOS SDP (13). Then the polynomial z\{x) := 1 + a'''(x)u* vanishes 
at the r+ points of the support of and the polynomial z*_{x) ■.= 1 — a'''(x)u* vanishes at 
the r_ points of the support of p*_. 

Proof: Let us denote by {x^}fc=i,...,r+ C X the points of the support of the optimal measure 
/i^, computed from the moments y'f solving optimally moment SDP (11). By complemen¬ 
tarity of the solutions of primal-dual SDP (11) and (13), it holds {z^, p*_) = 0 and hence 
{z\, 5^k^) = z*j_{x^) = 0 for each fc = 1,..., r+. The proof is similar for z*_ and p*_. □ 

4 Discussion 

We would like to point out that the developments in this paper were inspired by a previous 
work on optimal control for linear systems formulated as a primal LP (8) on measures and a 
dual LP on continuous functions (9), and solved numerically with primal-dual moment-SOS 
SDP hierarchies [CAHL13, CAHL14]. Formulating optimal control problems as moment 
problems was a classical research topic in the 1960s, where optimal control laws were sought 
in measures spaces (completions of Lebesgue spaces) to allow for oscillations and concentra¬ 
tions, see e.g. [Kra68] or the overview in [Fat99, Section III]. In the case of linear optimal 
control of an ordinary differential equation of order n, it was proved in [Neu64] that there is 
always an u-atomic optimal measure solving problem (8). 

In practice. Theorem 1 should be used as follows: 

1 Let k = max{d, kx}- 
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2 Solve SDP problem (11) and its dual (13) with a primal-dual algorithm. 

3 If the rank condition (12) of Theorem 1 is satished, then extract the measure from the 
solution of (11) and the polynomial certihcate from the solution of (13). Otherwise, 
let fc = /c -|- 1 , and go to 1 . 

We conjecture that if the data a, b in problem ( 8 ) are generic, then there is a hnite value of 
k for which the rank condition of Theorem 1 is satished. The rationale behind this assertion 
follows from a result by Nie [Niel4] on generic hnite convergence for the moment-SOS SDP 
hierarchy for polynomial optimization over compact basic semi-algebraic sets. Translated in 
the present context for a hxed family of data a, results in [Niel4] yield that there is a set of 
polynomials {hi,..., h^} C M[u], such that, given a feasible solution u of (9), if ^^(u) ^ 0 
for all £ = 1 ,..., L, then indeed 


nx 

1 + a^(x)u = Po(x) + gj{x), x e M*", 

i=i 

and 

nx 

1 - a^(x)u = pI{x) + '^P^ix) gj{x), x e M", 

1=1 

for some SOS polynomials , k = 1,2 and j = 1 ,..., nx- So if the optimal solution u* of (9) 
satishes ^^(u*) 7 ^ 0, £ = 1 ,..., L, then d* = dl for some index k (i.e. hnite convergence takes 
place). Similarly, by [Niel3] the rank-condition (12) of Theorem 1 also holds generically for 
polynomial optimization (which however is a context diherent from the present context). Put 
diherently, hnite convergence would not hold only if every optimal solution u of (9) would 
be a zero of some polynomial of the family {hi,..., h^} C M[u]. But so far we have not 
proved that at least one optimal solution u* of (9) is not a zero of some of the polynomials 
hi, at least for generic b. 

Of course, hnite convergence occurs for trigonometric polynomials on X = [0,27r], which 
follows from the Fejer-Riesz theorem and this was exploited in the landmark paper [CFG14]. 
Similarly, but apparently not so well-known, the Fejer-Riesz theorem also holds in dimension 
n = 2. Indeed it follows from Corollary 3.4 in [Sch06] that every non-negative bivariate 
trigonometric polynomial can be written as a sum of squares of trigonometric polynomials^. 
So again for trigonometric polynomials on X = [0,27r]^, hnite convergence of the hierarchy 
(13) takes place, i.e., dl = d*. Note however that in contrast to the one-dimensional case, 
there is no explicit upper bound on the degrees of the sum of squares which are required, so 
that even in the two-dimensional Fourier case we do not have an a priori estimates on the 
smallest value of k for which dl = d* an for which we can guarantee that the rank condition 
of Theorem 1 is satished. 

Generally speaking, even if our genericity conjecture is true, we do not have a priori estimates 
on the smallest value of k for which Theorem 1 holds. As mentioned above, this also true 
even in the two-dimensional case on [ 0 , 27r]^ where hnite convergence is guaranteed in all 
cases. 

^We are grateful to Markus Schweighofer for providing this reference. 
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5 Appendix 

We first recall some standard results of convex analysis. 


Lemma 6 ([FK94, Corollary 1.1.3]) Let C <Z ML he a closed convex cone with dual C* = 
{y : {x,y) >0, Vx G C}. Then int C* 7 ^ 0 (7 fl (—C) = {0}. 

Lemma 7 ([FK94, Corollary 1.1.6]) Let C C M" he a closed convex cone whose dual C* has 
nonempty interior. Then for all y G int C*, the set {x E C : {x,y) < 1} is compact. 

Lemma 8 The convex cone ^fc(gx) is closed. 


Proof: Let Sf be the convex cone of real symmetric matrices of size n that are positive 
semidehnite. Let be the set of n-dimensional integer vectors a such that oti <k and 
let Vk{x) := (x")aeN" be a vector of monomials of degree up to k. Next let Vk{x) Vk{x)'^ = 
^oa and 

Vk-v^ (x) Vk-vj (x)'^ gj{x) = ^ j = l,...,nx, 

for some appropriate real symmetric matrices Aja- 

Consider a sequence {qt)teN C ^fc(gx) such that qt ^ q E M[x]2A: as t —)■ cxd. That is 

nx 

qtix) = Potix) + ^Pjtix) gj{x), Vx G M”, 
i=i 

for some Pjt E T,[x]k-vj, j = 0,..., nx, for all t eN. More precisely, coefficient-wise 

nx 

qta = {Qot, Apa) + Aja), Vo G (14) 

i=i 

for some appropriate matrices Qjt E Let y = (|/Q,)Q,gNn^ be the moments ya '■= x°‘dx 

of the measure uniformly supported on X. Observe that since X has nonempty interior, 


p{x)dx >0 VO 7 ^pGS[x]fc, 


' X 


and 


' X 


p{x) gj{x)dx > D VO 7 ^ p G S[x]fc_^., j = l,...,nx- 


Put differently Mk{y) >- 0 and Mk-vj{gj y) >- 0, j = I,.. .,nx- 

The convergence qt ^ q implies {qt,y) —t {q,y) as t —)■ cxd. Hence there is some g such that 
V ^ {qt,y) for all t eN. This in turn implies 

nx 

V > {qt,y) = {pot,y)+ ^{Pjtgj,y) 


i=i 


nx 


{Qot, Mk{y)) + ^{Qjt, Mk-y.{gj y)). 
i=i 


(15) 
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Therefore 


sup (Qot, Mkiy)) < r], sup{Qjt, Mk-v^ {gj y)) < g, j = 1,..., nx- 


As 0 -< Mk{y) G and 0 -< Mk-v^igjy) G j = one may 

invoke Lemma 7 and conclude that the sequences (QoOieN C and {Qjt)t&n C 
are norm-bounded. Therefore there is a subsequence {te)£^^ and matrices Qq G and 
Qj ^ S’l j = 1,..., nx, such that 


Qoti — t Qo, Qjti Qj, j — 1, • • •, u-x 


as £ —)■ oo. Taking the limit for the subsequences {qtia)ieN uud {Qjtjien in (14) yields 
coefficient-wise 


*?« {Qo, ^Oa) T ^ ^ {Qjy ^ja )) 



i=i 

which proves that q G ^k{Sx), the desired result. 


□ 
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